Ground-state behavior of the 3d ±J random-bond Ising model 
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Large numbers of ground states of the three-dimensional 
± J random-bond Ising model are calculated for sizes up to 14 3 
using a combination of a genetic algorithm and Cluster-Exact 
Approximation. Several quantities are calculated as function 
of the concentration p of the antiferromagnetic bonds. The 
critical concentration where the ferromagnetic order disap- 
pears is determined using the Binder cumulant of the magne- 
tization. A value of p c = 0.222 ± 0.005 is obtained. From the 
finite-size behavior of the Binder cumulant and the magneti- 
zation critical exponents v = 1.1 ± 0.3 and /3 = 0.2 ± 0.1 are 
calculated. 

The behavior of the distribution of overlaps P(q) is used to 
investigate how the spin-glass phase evolves with increasing 
concentration p. The spin-glass order is characterized by a 
broad distribution of overlaps which extends down to q = 0. 

Keywords (PACS-codes): Spin glasses and other 
random models (75.10.Nr), Numerical simulation studies 
(75.40.Mg), General mathematical systems (02.10.Jf). 

Introduction In this work systems of N spins o~i = 
±1, described by the Hamiltonian 
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are investigated. The spins are placed on a three- 
dimensional (d=3) cubic lattice of linear size L with 
periodic boundary conditions in all directions. Systems 
with quenched disorder of the nearest-neighbor interac- 
tions (bonds) are investigated. Their possible values are 
Jij = ±1. The concentration of the antiferromagnetic 
(AF) bonds (Jy = —1) is denoted withp, all other (1— p) 
interactions are ferromagnetic. A constrained disorder is 
used, so that the fraction of the antiferromagnetic bonds 
is exactly p for all realizations of the disorder. 

The model shows a complex behavior for low temper- 
atures. For large concentrations of the ferromagnetic 
bonds it is energetically favorable for two interacting 
spins to have the same orientation. So the system shows 
ferromagnetic order, which means that most of the spins 
have the same value. For large concentrations p the sys- 
tem is antiferromagnetically ordered: the system can be 
divided into to penetrating sublattices and each sublat- 
tice has ferromagnetic order, but the sign of the ordering 
of the two sublattices is different. For intermediate con- 
centrations of the antiferromagnetic bonds neither fer- 
romagnetic nor antiferromagnetic order exists. The sys- 
tem is called a spin glass |ij . For finite-dimensional spin 
glasses no final agreement about their behavior exists. 



Recent results from simulations [gH of small systems 
with p = 0.5 up to size 16 3 indicate that the three- 
dimensional spin glass has a complex behavior: the (free) 
energy landscape has many stable configurations which 
differ strongly from each other. Whether the onset of this 
spin-glass behavior takes place at the same concentration 
where the ferromagnetic order disappears is unclear. 

Here the ground-state, i.e. zero temperature (T = 0), 
behavior of the model as function of the concentration p 
is investigated. Since the phase diagram of the system is 
symmetrical to p — 0.5 if one identifies the ferromagnetic 
with the antiferromagnetic regime, only values p < 0.5 
are used. Of special interest is the critical concentration 
p c where the ferromagnetic order disappears. 

In the past the ±J random-bond Ising model has 
been studied using Monte-Carlo simulations Q, zero- 
temperature series expansions 0, high-temperature se- 
ries expansions , Monte Carlo renormalization-group 
calculations Q and renormalization-group theory ||. 
All of these results are qualitatively consistent with the 
phase-diagram described above, but no agreement is 
found on the value of the critical concentration p c or 
its temperature dependence. In all of these publications 
the detailed structure of the states was not investigated, 
which can be done for example by calculating the distri- 
bution of overlaps. 

A study of this model using true ground states has not 
been performed before. Only for special two-dimensional 
systems, where exact ground states can be calculated ef- 
ficiently, some results [[To] |l2| are known. There the crit- 
ical concentration of a square lattice is estimated to be 
pf = 0.10. 

The behavior of the ± J random bond Ising model is 
determined by the occurrence of frustration ]l3{ . The 
simplest example of a frustrated system is a triple of 
spins where all pairs are connected by antiferromagnetic 
bonds. It is not possible to find a spin-configuration were 
all bonds contribute with a negative value to the energy. 
One says that it is not possible to satisfy all bonds. In 
general a system is frustrated, if closed loops of bonds 
exists, where the product of these bond-values is neg- 
ative. For square and cubic systems the smallest closed 
loops consist of four bonds. They are called (elementary) 
plaquettes. 

The presence of frustration makes the calculation of ex- 
act ground states of such systems computationally hard. 
Only for the special case of the two-dimensional system 
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with periodic boundary conditions in no more than one 
direction and without external field a polynomial-time al- 
gorithm is known |Q . For the general case the simplest 
method works by enumerating all 2 N possible states and 
has obviously an exponential time complexity. Even a 
system size of 4 3 is too large. The basic idea of an al- 
gorithm called Branch- and- Bound [ fl5| is to exclude the 
branches of the search-tree, where no low-lying states can 
be found, so that the complete low-energy landscape of 
systems of size 4 3 can be calculated Jl6| . A more sophis- 
ticated method called Branch- and- Cut |lj]|l8| works by 
rewriting the problem as a linear optimization problem 
with an additional set of inequalities which must hold for 
the solution. Since not all inequalities are known a priori 
the method iteratively solves the linear problem, looks 
for inequalities which are violated, and adds them to the 
set until the solution is found. Since the number of in- 
equalities grows exponentially with the system size the 
same holds for the computation time of the algorithm. 
With Branch-and-Cut anyway small systems up to 8 3 
are feasible. The method used here is able to calculate 
true ground states up to size 14 3 . 

By studying ground states one does not encounter er- 
godicity problems or critical slowing down like in Monte- 
Carlo simulations. Since it is possible to compare with 
exact results from Branch-and-Cut calculations no un- 
controlled approximations are used. The only uncer- 
tainty comes from the fact that one is restricted to rela- 
tively small systems. 

In the next section the algorithm is explained. Then 
all results are presented. In the last section a conclusion 
is driven. 



ALGORITHM 

The algorithm for the calculation of the ground states 
bases on a special genetic algorithm 19 2(]] and on 
Cluster-Exact Approximation (CEA) |2lf| which is a so- 
phisticated optimization method. Next a short sketch of 
these algorithms is given. 

The genetic algorithm starts with an initial popula- 
tion of Mi randomly initialized spin configurations (= 
individuals) , which are linearly arranged in a ring. Then 
n/j x Mi times two neighbors from the population are 
taken (called parents) and two offsprings are created us- 
ing a triadic crossover: a mask is used which is a third 
randomly chosen (usually distant) member of the pop- 
ulation with a fraction of 0.1 of its spins reversed. In 
a first step the offsprings are created as copies of the 
parents. Then those spins are selected, where the orien- 
tations of the first parent and the mask agree The 
values of these spins are swapped between the two off- 
springs. Then a mutation with a rate of p m is applied to 
each offspring, i.e. a fraction p m of the spins is reversed. 

Next for both offsprings the energy is reduced by ap- 
plying CEA: The method constructs iteratively and ran- 



domly a non-frustrated cluster of spins. Spins adjacent 
to many unsatisfied bonds are more likely to be added to 
the cluster. During the construction of the cluster a local 
gauge-transformation of the spin variables is applied so 
that all interactions between cluster spins become ferro- 
magnetic. 




FIG. 1. Example of the Cluster- Exact Approximation 
method. A part of a spin glass is shown. The circles repre- 
sent lattice sites/spins. Straight lines represent ferromagnetic 
bonds the jagged lines antiferromagnetic interactions. The 
top part shows the initial situation. The construction starts 
with the spin at the center. The bottom part displays the fi- 
nal stage. The spins which belong to the cluster carry a plus 
or minus sign which indicates how each spin is transformed, 
so that only ferromagnetic interactions remain inside the clus- 
ter. All other spins cannot be added to the cluster because it 
is not possible to multiply them by ±1 to make all adjacent 
bonds positive. Please note that many other combinations of 
spins can be used to build a cluster without frustration. 

Fig. [j] shows an example of how the construction of 
the cluster works using a small spin-glass system. For 3d 
± J spin glasses each cluster contains typically 58 percent 
of all spins. The non-cluster spins act like local magnetic 
fields on the cluster spins, so the ground state of the clus- 
ter is not trivial. Since the cluster has only ferromagnetic 
interactions, an energetic minimum state for its spins can 
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be calculated in polynomial time by using graph theo- 
retical methods [|23H25f : an equivalent network is con- 
structed |26|, the maximum flow is calculated ]27| , |2l| p1 
and the spins of the cluster are set to their orientations 
leading to a minimum in energy. This minimization step 
is performed n m j n times for each offspring. 

Afterwards each offspring is compared with one of its 
parents. The pairs are chosen in the way that the sum 
of the phenotypic differences between them is minimal. 
The phenotypic difference is defined here as the number 
of spins where the two configurations differ. Each parent 
is replaced if its energy is not lower (i.e. not better) 
than the corresponding offspring. After this whole step is 
done nfj x Mi times, the population is halved: From each 
pair of neighbors the configuration which has the higher 
energy is eliminated. If more than 4 individuals remain 
the process is continued otherwise it is stopped and the 
best individual is taken as result of the calculation. 

The representation in fig. || summarizes the algorithm. 

The whole algorithm is performed hr times and all 
configurations which exhibit the lowest energy are stored, 
resulting in fig statistically independent ground-state 
configurations. 

This algorithm was already applied to examine the 
ground state structure of 3d spin glasses j| . 



1 Implementation details: We used Tarjan's wave algorithm 
together with the heuristic speed-ups of Traff. In the con- 
struction of the level graph we allowed not only edges (v,w) 
with level(u)) = level (v)+l, but also all edges (v,t) where t is 
the sink. For this measure, we observed an additional speed- 
up of roughly factor 2 for the systems we calculated. 



algorithm genetic CEA({Jy}, Mi, n R , p m , n min ) 
begin 

create Mi configurations randomly 

while (Mi > 4) do 

begin 

for i = 1 to tir x Mi do 
begin 

select two neighbors 

create two offsprings using triadic crossover 
do mutations with rate p m 
for both offsprings do 
begin 

for j = 1 to n m i n do 

begin 

construct unfrustrated cluster of spins 
construct equivalent network 
calculate maximum flow 
construct minimum cut 
set new orientations of cluster spins 
end 

if offspring is not worse than related parent 
then 

replace parent with offspring 

end 
end 

half population; M{ = Mi/2 
end 

return one configuration with lowest energy 
end 

FIG. 2. Genetic Cluster-exact Approximation. 



RESULTS 

We used the simulation parameters determined in for- 
mer calculations for p — 0.5: For each system size 
many different combinations of the simulation parame- 
ters mi, riR, n m i n ,p m were tried for some sample systems. 
The final parameters where determined in a way, that by 
using four times the numerical effort no reduction in en- 
ergy was obtained. Here p m = 0.2 and tlr = 10 were 
used for all system sizes. 

For smaller concentrations p the ground states are eas- 
ier to find, because the number of frustrated plaquettes 
is smaller. But it was not possible to reduce the compu- 
tational effort substantially in order to get still ground 
states. So we used the parameters of p = 0.5 for all con- 
centrations p. Table 1 summarizes the parameters. Also 
the typical computer time r per ground state computa- 
tion on a 80 MHz PPC601 is given. Using these param- 
eters on average ng > 8 ground states were obtained for 
every system size L using n R = 10 runs per realization. 
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Tab 1. Simulation parameters: L — system size, Mi 
= initial size of population, tir — average number of off- 
springs per configuration, n m i n = number of CEA mini- 
mization steps per offspring, r = average computer time 
per ground state on a 80MHz PPC601. 

We compared our results for 180 sample systems of L = 
6 with exact ground states which were obtained using 
a Branch-and-Cut program The genetic CEA 

algorithm found the true ground states for all systems! 
The same result was obtained for L = 4 as well.^J A 
more detailed analysis is presented in (29). So we can be 
sure that genetic CEA and our method of choosing the 
parameters lead to true ground states or at least to states 
very close to true ground states. 

We performed ground state calculations for p G 
[0.1,0.4] for lattice sizes L = 3,4,5,6,8,10,14. The 
number of independent realizations of the bond-disorder 
ranged from N L = 1000 for L = 14 to 30000 for smaller 
systems. Most of the systems have AF-bond concentra- 
tions around p = 0.22 



The ground state energy e as function of system size 
for different system sizes is shown in Fig |. To keep the 
figure clear only the sizes L — 3,4,6,14 are presented. 
For small concentrations p the ground state is mainly 
ferromagnetic. It follows that all AF bonds are not sat- 
isfied, so the ground state energy increases linearly like 
e(p) « —3 + 6p with p. For larger concentrations the 
ground state energy approaches the p — 0.5 limit, be- 
cause the spins can arrange so that not all AF-bonds are 
broken. With increasing L the ground state energy de- 
creases, because the periodic boundary conditions impose 
less constraints on the system. For L — > oo and p = 0.5 
a ground state energy of eoo(0.5) = —1.7876(3) is found 
m§. 
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FIG. 3. Average ground state energy e per spin as function 
of AF-bond concentration p for system sizes L — 3, 4, 6, 14. 
Lines are guides for the eyes only. 



2 For L > 6 the Branch-and-Cut program needs to much 
computer time because of the exponential time complexity. 



FIG. 4. Binder cumulant of magnetization as function of 
AF-bond concentration p for system size L = 3,5, 8, 10, 14. 
Two typical error bars are given. Lines are guides for the 
eyes only. 

From fig. H is clear that the energy as function of con- 
centration is not well suited for determining the critical 
concentration p c where the ferromagnetic behavior dis- 
appears. For this purpose the Binder cumulant []30| , j3l]| 

*' x » s K 3 -(w) (2) 

for the magnetization M = -h Oi is used. The av- 
erage (. . .) denotes both average over different ground 
states of a realization and over the disorder. In fig. |J the 
Binder cumulant is shown for L = 3, 5, 8, 10, 14. L=4, 6 
are omitted in this figure to keep it clear. For the same 
reason only typical error bars for two sample points are 
shown. All curves intersect at p c = 0.222 ± 0.002. Only 
L = 4 (not shown) is little worse, because it meets the 
others in the interval p £ [0.217,0.222]. So we conclude 
that the critical concentration for the ferromagnetic or- 
der is p c = 0.222(5). 
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The value for p c is comparable to results from high- 
temperature series expansions: p c = 0.19(2) p c as 
0.25 0, from Monte-Carlo renormalization-group results 
p c = 0.233(4) m and from Monte-Carlo simulations: 
p c « 0.24 0. Our value is much larger than a result from 
a zero-temperature expansion: p c — 0.12 — 0.13 (^] and 
much lower than a recent result from a renormalization- 
group study: p c « 0.37 ||. Since the intersection of the 
curves of the Binder cumulant is very sharp, we believe 
that our result is very reliable, although the systems in- 
vestigated here are rather small. 




-0.3 -0.2 
(p-0.222) ' 



FIG. 5. Scaled plot of Binder cumulant. Line is guide for 
the eyes only. 

For the Binder cumulant the following finite-size scal- 
ing relation is assumed pW 



q(p,L) = ~g(L 1 /»(p-pc)) 



(3) 



By plotting g(p, L) against (p—p c ) with correct pa- 
rameter v the datapoints for different system sizes should 
collapse onto a single curve. The best results were ob- 
tained for pc = 0.222 and l/v = 0.9. In fig. | the 
resulting scaling plot is shown. It is possible to change 
the value of v in a wide range without large effects on the 
scaling plot. So we estimate v = 1.1(3). This is consis- 
tent with v = 1.7(3) which was found using Monte-Carlo 
simulations of spin glasses (p — 0.5) at finite temperature 

H- 

The average magnetization m = (M) has the standard 
finite-size scaling form §33] 



m(p, L) = L- f3/y m(L 1/u {p - p c )) 



(4) 



By plotting L^/ V m{p, L) against L x l h '(p — p c ) with cor- 
rect parameters (3, v the datapoints for different system 
sizes should collapse onto a single curve. The best re- 
sult was obtained using \ jv = 0.9 and (5/v = 0.19. It 



is shown in fig. || for L = 3,5, 8, 10, 14. From variations 
of the value (3/v we estimate the value of the exponent 
= 0.2(1). 




o.o 

(p-0.222) * L" 

FIG. 6. Scaled plot of magnetization. The inset show the 
raw data for L = 3,5, 14. Lines are guides for the eyes only. 

To characterize the spin-glass behavior the overlap q is 
used. It compares two different states {erf}, {uf } of the 
same realization of the random bonds 
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For ferromagnetic/antiferromagnetic order two inde- 
pendently calculated ground states are identical or re- 
lated by a global flip of all spins, i.e. q = ±1. For 
spin-glass order many different ground states exist |s| , so 
q can take also intermediate values g G [—1,1]. Since the 
system has no external field, each state is equivalent to 
the state where all spins are reversed, so only the absolute 
value \q\ is considered here. 
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FIG. 7. Average overlap (\q\) value as function of AF-bond 
concentration p for L = 3,4, 8, 14. With increasing concentra- 
tion more and more spins belong to clusters which contribute 
to the degeneracy of the ground state, so the average overlap 
value decreases. Where the onset of the spin glass behavior 
is located exactly cannot be seen from this figure. Lines are 
guides for the eyes only. 

In fig. [t] the average value of the overlap (\q\) = (\q a>3 \} 
is shown for the lattice sizes L — 3,4, 8, 14. The decrease 
of with increasing concentration p is clearly visible. 
But this quantity has much larger fluctuations than the 
magnetization, so it is difficult to use this data for fur- 
ther analysis. Also the value is not equal to 1.0 for a large 
interval, so it is not possible to extract at which concen- 
tration the onset of the spin-glass behavior is located. 




FIG. 8. Distribution of overlaps for AF-bond con- 

centrations about p « 0.18 for L = 3, 4, 8, 14. A finite frac- 
tion of spins is contained in small clusters which can take two 
orientations in the ground state. For L — > oo a delta-function 
is obtained similar to the distribution found for the ground 
states of random-field Ising systems. Lines are guides for the 
eyes only. 

More information can be obtained, if one calculates 
not only the average of \q\, but its distribution [B4| 



P(\q\) = (S(\q\ - |<r|)> 



(G) 



which describes the ground-state structure. In fig. g 
the distributions for sizes L = 3,4, 8, 14 at p w 0.18 are 
presented]^]. With increasing size the distributions be- 
comes narrower. The reason is, that due to the discrete 



bond distribution Jy = ±1 there are always some small 
clusters of spins, which can take two orientations in the 
ground state. With increasing system size L these ef- 
fects cancel out and P (\q\) converges to a delta- function 
S(q — qf ree ), where g/ ree is just the fraction of spins con- 
tained in such free clusters. This ground-state structure 
is similar to that of random-field Ising systems [p5| . 




FIG. 9. Distribution of overlaps for AF-bond con- 

centrations about p » 0.23 for L — 3, 4, 8, 14. The distribu- 
tion is broad and extends to q — for all sizes L, indicating 
a spin-glass behavior. Lines are guides for the eyes only. 

In fig. ^ the distributions of overlaps is displayed for a 
concentration slightly larger than p c . Here the behavior 
is completely different. The distributions extend over a 
large interval down to q = 0. With increasing lattice 
size L only the shape of the large-q part changes a little 
bit while for small values no systematic modification is 
visible. The peak at large q- values, which raises with 
increasing concentration p, results from small clusters of 
spins which can take two orientations in the ground state. 
This is the same as for smaller concentrations of the AF- 
bonds. But the large extent down to q = cannot be 
explained in this way. It shows that spin-glass ground 
states have a very rich structure, similar to the behavior 
found for the SK- model fl34j, where each spin interacts 
with every other spin (for a detailed discussion of the 
ground-state structure see p|,p6[). 

To investigate this behavior quantitatively the vari- 
ance cr 2 (|<z|) of the distributions as function of AF-bond 
concentration p and system size L is calculated 



q\) = ((\q ap \-(\q\)r) 



(7) 



3 Since all realizations of a given size L have exactly the same 
number of ferromagnetic bonds, for each system size only a each lattice size the value of p is chosen which is nearest to 
finite number of different concentrations is possible. Here for 0.18 
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FIG. 10. Average variance o" 2 (|g|) of overlap distribution 
as function of AF-bond concentration p. For small concen- 
trations the width of the distribution shrinks to zero with 
increasing size L. For larger values the spin-glass phase is 
characterized by broad distributions of overlaps. Lines are 
guides for the eyes only. 

In fig. |l0| the result for L = 3,4,8,14 is shown. For 
small concentration p the width of the distributions is 
small and shrinks with increasing size L. For larger con- 
centrations the width increases and remains nonzero even 
for larger lattice sizes. In || it was shown, that the spin 
glass ground state with p = 0.5 is likely to have a broad 
distribution even for L — > oo. For larger sizes the statis- 
tics is too bad to extract for example a critical concentra- 

(2) 

tion p c by a finite-size scaling analysis similar to that 
presented above. 

The onset of the spin glass behavior can even better 
be observed by calculating the fraction X qo of the distri- 
bution of overlaps below a fixed value qo : 

X qo = / P(\q\)dq (8) 
JO 




FIG. 11. Average fraction Xo.5 of overlap distribution be- 
low q = 0.5 as function of AF-bond concentration p for sizes 
L — 3, 4, 8, 14. In the spin-glass phase the distribution of 
overlaps extends to q — for all lattice sizes. Lines are guides 
for the eyes only. 

In fig. [ll] the value of A"o.5 is shown as function of p 
for the lattice sizes L = 3, 4, 8, 14. For the limiting case 
p = 0.5 the value of X 5 converges to a nonzero value 
for L -> 00 [0. 

CONCLUSION 

Using a combination of a genetic algorithm and 
Cluster-Exact Approximation ground states of the three- 
dimensional ±J random-bond Ising model were calcu- 
lated for different concentrations of the antiferromagnetic 
bonds. A former comparison with exact ground states 
calculated using a Branch-and Cut program shows that 
genetic CEA is able to calculate true ground states. 

For small concentrations the ground state is mainly 
ferromagnetic. The critical concentration where the fer- 
romagnetic order disappears was determined using the 
Binder-cumulant g(p, L) of the magnetization: 

p c = 0.222 ± 0.005 (9) 

This is the first time p c is calculated by direct investi- 
gations of ground states. 

Using a finite-size scaling analysis of the magnetiza- 
tion and the Binder-cumulant critical exponents were ob- 
tained: 

v = 1.1 ± 0.3, /3 = 0.2 ± 0.1 (10) 

These values are consistent with results from Monte- 
Carlo simulations of the p = 0.5 case at finite temper- 
ature, where the same analysis presented here was per- 
formed for the spin-glass order parameter q instead of the 
magnetization m. 



7 



The onset of the spin-glass behavior with increasing 
lattice size was investigated the first time by calculat- 
ing the distributions of overlaps as function of p. The 
spin glass phase is characterized by a broad distribution 
of overlaps which extends down to q — and does not 
change substantially with increasing system size. For this 
quantity the bad statistics for larger values of the con- 
centration p > 0.23 does not allow to determine a second 

(2) 

critical concentration p c , so it is yet not possible to 
check whether the ferromagnetic order disappears at the 
same concentration where the spin-glass phase appears. 
Here much more data are needed. 
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